Introduction

Bulk RNA-seq analysis of cerebral organoid samples in conditions co-cultured (with/without microglia) and stimulation (ctr/LPS/IFNy).

CO vs COiMG

Organoids co-cultured with vs without microglia.

Variance analysis

PCAplot(pca_cor3_mod1, "COiMg", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now

ggplot(pca.df_upgraded,aes_string(x='COiMG',y='Phagocytosis', fill='COiMG')) +
  stat_compare_means(method = "t.test", label.y = 10, label.x = 1.3) +
  geom_boxplot() +
  geom_boxplot(width=0.1, fill="white") +
  scale_fill_brewer(palette="Dark2") +
  labs(title="Phagocytosis differential PC expression",x="COiMG", y = "PC1 Phagocytosis") +
  theme_classic()

Heatmap

#for control only
ctr.meta <- meta3[meta3$condition %in% 'ctrl',]
ctr.meta$COiMg[ctr.meta$COiMg %in% 'yes'] <- 'COiMg'
ctr.meta$COiMg[ctr.meta$COiMg %in% 'no'] <- 'CO'

#horizontal
ra <- rowAnnotation(
  COiMG = ctr.meta$COiMg[ctr.meta$condition %in% 'ctrl'],
  col = list(
    COiMG = c("COiMg"='tomato','CO'='lightblue')),
  show_annotation_name = T,
  show_legend = T)


ctr <- t(batch.rem3_mod1[,meta3$condition %in% 'ctrl'])
hm_counts <- scale(ctr[,colnames(ctr) %in% gene_panels_subset$Microglia])
ComplexHeatmap::Heatmap(hm_counts,
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = T,
                        show_column_dend = T,
                        show_row_names = F,
                        show_column_names = F,
                        left_annotation = ra,
                        bottom_annotation = HeatmapAnnotation(
                        text = anno_text(colnames(hm_counts), rot = 300, just = "left",gp=gpar(fontsize=8))),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
                        name = "Z-score expression")

DEA

Volcano

volcano_plot <- function(res, title = NULL, subtitle = NULL, annotate_by = NULL, type ='ALS', ymax1 = 7.5, ymax2 = 8, xmax1 = -4.5, xmax2 = 4.5){
  res <- 
    mutate(res,
           sig = case_when(
             padj >= 0.05 ~ "non_sig",
             padj < 0.05 & abs(log2FoldChange) < 1 ~ "sig",
             padj < 0.05 & abs(log2FoldChange) >= 1 ~ "sig - strong"
           )) %>%
    mutate(direction = ifelse(log2FoldChange > 0, "up", "down")) %>%
    mutate(log2FoldChange = case_when(
      log2FoldChange > 10 ~ Inf,
      log2FoldChange < -10 ~ -Inf,
      TRUE ~ log2FoldChange
    )) %>%
    mutate(class = paste(sig, direction))
  if( type == "ALS"){
    xpos <- 0.5
    ymax <- ymax1
    xlim <- c(xmax1,xmax2)
  }else{
    xpos <- 0.025
    ymax <- ymax2
    xlim <- c(-0.042, 0.042)
  }
  de_tally <- group_by(res, sig, direction, class) %>% tally() %>%
    filter(sig != "non_sig") %>%
    mutate( position = ifelse(sig == "sig", xpos, 2) ) %>%
    mutate(position = ifelse( direction == "down", -1 * position, position)) %>%
    mutate(n = formatC(n, format="f", big.mark=",", digits=0))
  plot <- res %>%
    mutate( pvalue = ifelse( pvalue < 1e-90, Inf, pvalue)) %>% #threshold at 1e16
    ggplot(aes(x = log2FoldChange, y = -log10(pvalue))) + 
    #geom_point(aes(colour = class ), size = 0.5) +
    rasterise(geom_point(aes(colour = class ), size = 0.5), dpi = 300) +
    scale_colour_manual(values = c("non_sig up" = "gray", 
                                   "non_sig down" = "gray",
                                   "sig up" = "#EB7F56", 
                                   "sig - strong up" = "#B61927",
                                   "sig down" = "#4F8FC4",
                                   "sig - strong down" = "dodgerblue4"
    )) +
    theme_bw() +
    labs(y = expression(-log[10]~P~value), x = expression(log[2]~"(fold change)"), title = title, subtitle = subtitle) +
    guides(colour = FALSE) +
    scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) +
    theme_bw() +
    theme(
      plot.title = element_text(hjust = 0.5),
      plot.subtitle = element_text(hjust = 0.5),
      panel.border = element_blank(),
      axis.ticks = element_line(colour = "black")
    ) +
    geom_text(fontface = "bold", data = de_tally, aes(x = position, y = ymax - 0.5, label = n, colour = class), size = 4) +
    scale_x_continuous(limits = xlim)
  if(!is.null(annotate_by)){
    plot <- plot + 
      ggrepel::geom_text_repel(
        fontface = "italic",
        data = filter(res, symbol %in% annotate_by), 
        aes(x = log2FoldChange, y = -log10(pvalue), label = symbol), 
        min.segment.length = unit(0, "lines"),
        size = 2.3) +
      geom_point(
        data = filter(res, symbol %in% annotate_by), size = 0.8, colour = "black"
      ) +
      geom_point(aes(colour = class ),
                 data = filter(res, symbol %in% annotate_by), size = 0.6
      )
  }
  return(plot)
}
res_COiMG$symbol <- rownames(res_COiMG)
volcano_plot(data.frame(res_COiMG), title = 'CO vs COiMg',
             annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)

Top 50 DEGs heatmap

top15.coimg <- c(rownames(results$COiMG[order(results$COiMG$log2FoldChange),][results$COiMG[order(results$COiMG$log2FoldChange),]$padj < 0.05,][1:25,]),
rownames(results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,]))

ha <- HeatmapAnnotation(
                        COiMG = meta3$COiMg,
                        col = list(
                                    COiMG = c("yes"='tomato','no'='lightblue')))
ComplexHeatmap::Heatmap(t(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top15.coimg,]))),
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        show_column_names = F,
                        top_annotation = ha,
                        row_names_gp = gpar(fontsize=8),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(30),
                        name = "Z-score expression",
)

Which genes of the top DEGs are cilia-associated (CBC)?

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)

organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
                               'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]


knitr::kable(cbind('DEG'=top15.coimg,'Cilia-associated'=top15.coimg %in% organoid$CBC))
DEG Cilia-associated
CD164L2 FALSE
STOML3 TRUE
MMP10 FALSE
GIPC2 FALSE
LRRC18 FALSE
FAM166C FALSE
LRRC71 FALSE
IL5RA FALSE
TTC22 FALSE
CD302 FALSE
C20orf85 FALSE
CLDN2 FALSE
MFRP FALSE
TACR3 FALSE
CIBAR2 FALSE
C11orf97 FALSE
F5 FALSE
SAG FALSE
ADGB FALSE
TEKT4 FALSE
VWA5B1 FALSE
DEUP1 FALSE
DYNLT5 FALSE
ARSI FALSE
FAM81B TRUE
APOC2 FALSE
GOLGA6D FALSE
VAV1 FALSE
FCGR3A FALSE
C1QC FALSE
ITGB2 FALSE
CD86 FALSE
OLR1 FALSE
SPP1 FALSE
WDFY4 FALSE
ALOX5AP FALSE
TYROBP FALSE
ADAMDEC1 FALSE
BIN2 FALSE
CD53 FALSE
MS4A7 FALSE
SRGN FALSE
CCRL2 FALSE
CD52 FALSE
FCER1G FALSE
LSP1 FALSE
PECAM1 FALSE
CD48 FALSE
DCSTAMP FALSE
CHI3L1 FALSE

GO

dotplot(all_res$COiMG$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")

dotplot(all_res$COiMG$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")

Enrichment

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)

organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
                               'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]


#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
  DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))


for (x in colnames(organoid)){
  DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}

DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
                    "Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)

DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

For microglia pathways.

DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
  DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))


for (x in colnames(gene_panels_subset)){
  DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}

DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
                    "Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
cgenes <- readxl::read_xlsx('Gene lists for organoid paper.xlsx', col_names = T)
m <- data.frame(t(toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`))[toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`)) %in% rownames(bulkorg_counts3)]), check.names = F)
for(i in names(cgenes)[-1]){
  k <- na.omit(cgenes[i])
  s <- k[as.matrix(k)[,1] %in% rownames(bulkorg_counts3),]
  m <- rbind.fill(data.frame(m),data.frame(t(s)))
}
rownames(m) <- names(cgenes)
m <- t(m)
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
  ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))


for (x in colnames(m)){
  ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}

ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
                    "Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)

ComplexHeatmap::Heatmap(ndd.enrich.ress,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

IFNy stimulation

Let’s do this for IFNy stimulation

Variance analysis

PCAplot(pca_cor3_mod1, "condition", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now

ggplot(pca.df_upgraded,aes_string(x='Condition',y='`IFNy`', fill='Condition')) +
  #stat_compare_means(method = "t.test", label.y = 20, label.x = 2) +
  geom_boxplot() +
  geom_boxplot(width=0.1, fill="white") +
  scale_fill_brewer(palette="Dark2") +
  labs(title="Condition differential PC expression",x="Condition", y = "PC1 IFNy response") +
  theme_classic()

DEA

Volcano

res_IFNg$symbol <- rownames(res_IFNg)

volcano_plot(data.frame(res_IFNg), title = 'CTR vs IFNy',
             annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)

COiMg.degs <- all_res$COiMG$DEGs

Top 50 DEGs heatmap

top15.ifng <- c(rownames(results$IFNg[order(results$IFNg$log2FoldChange),][results$IFNg[order(results$IFNg$log2FoldChange),]$padj < 0.05,][1:25,]),
rownames(results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,]))

ha <- HeatmapAnnotation(Condition = meta3$condition,
                        COiMG = meta3$COiMg,
                        Line = meta3$Line,
                        col = list( Condition = c("ctrl" = "lightblue", "IFNg" = "tomato", "LPS" = "darkred"),
                                    COiMG = c("yes"='darkblue','no'='grey80'),
                                    Line = c('WTC11'='black','MSN38'='white')))
ComplexHeatmap::Heatmap(t(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top15.ifng,]))),
                        cluster_rows = T,
                        cluster_columns = T,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        show_column_names = F,
                        top_annotation = ha,
                        row_names_gp = gpar(fontsize=8),
                        col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(30),
                        name = "Z-score expression",
)

Which genes of the top DEGs are cilia-associated (CBC)?

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)

organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
                               'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]


knitr::kable(cbind('DEG'=top15.ifng,'Cilia-associated'=c(top15.ifng %in% organoid$CBC)))
DEG Cilia-associated
CROCC2 FALSE
F5 FALSE
STMND1 FALSE
MUSK FALSE
ARHGAP36 FALSE
PTPRC TRUE
IGF1 FALSE
KL FALSE
STOML3 TRUE
ATP13A5 FALSE
SLC39A12 TRUE
CYTL1 FALSE
C11orf97 FALSE
SERPINA3 FALSE
TMEM125 FALSE
FAM81B TRUE
C20orf85 FALSE
CAPSL TRUE
CIBAR2 FALSE
TRPM1 FALSE
C10orf105 FALSE
HMGCS2 FALSE
CFAP73 FALSE
PMCH FALSE
MMP10 FALSE
IDO1 FALSE
GBP5 FALSE
ETV7 FALSE
CXCL10 FALSE
SIGLEC7 FALSE
GBP1 FALSE
CXCL9 FALSE
APOL6 FALSE
APOL1 FALSE
PLAAT4 FALSE
GBP3 FALSE
SAMD9L FALSE
GBP2 FALSE
XAF1 FALSE
PSMB9 FALSE
CXCL11 FALSE
GBP6 FALSE
CD274 FALSE
OAS2 FALSE
DES FALSE
UBD FALSE
OR2I1P FALSE
OAS1 FALSE
NLRC5 FALSE
BATF2 FALSE

GO

dotplot(all_res$IFNg$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")

dotplot(all_res$IFNg$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")

Enrichment

For organoid cell-types.

setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)

organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
                               'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]


#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
  DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))


for (x in colnames(organoid)){
  DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}

DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
                    "Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)

DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

For NDD gene lists.

setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')

#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
  ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))


for (x in colnames(m)){
  ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}

ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
                    "Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)

ComplexHeatmap::Heatmap(ndd.enrich.ress,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)

DEGs$`Down-regulated`[DEGs$`Down-regulated` %in% as.character(na.omit(m[,12]))]
## [1] "CYP7B1" "GRIN2A"

For microglia pathways.

DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
  DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))


for (x in colnames(gene_panels_subset)){
  DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}

DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])

#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
                    "Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
                        cluster_rows = F,
                        cluster_columns = F,
                        show_row_dend = F,
                        show_column_dend = F,
                        show_row_names = T,
                        col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
                        name = "log2(q)",
)